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We have studied preheating of field perturbations in a 3-dimensional lattice including the effect 
of scalar metric perturbations, in two generic models of inflation: chaotic inflation with a quartic 
potential, and standard hybrid inflation. We have prepared the initial state for the classical evolution 
of the system with vanishing vector and tensor metric perturbations, consistent with the constraint 
equations, the energy and momentum constraints. The non-linear evolution inevitably generates 
■ vector and tensor modes, and this reflects on how well the constraint equations are fulfilled during 

the evolution. The induced preheating of the scalar metric perturbations is not large enough to 
backreact onto the fields, but it could affect the evolution of vector and tensor modes. This is the 
case in hybrid inflation for some values of the coupling g and the height of potential Vq^. For 
|. example with ~ 10 15 GeV, preheating of scalar perturbations is such that their source term 

in the evolution equation for tensor and vector fluctuations becomes comparable to that of the field 
00 ' anisotropic stress. 
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5-H ■ I. INTRODUCTION 

Cosmological observations, and in particular Cosmic Microwave Background (CMB) measurements [H, [H, [H, 0, [j| , 
are consistent with an early period of inflation, which among other things accounts for the inferred flatness of the 
Universe, and gives rise to the primordial curvature perturbation which would seed the large scale structure observed 
today. The value of the primordial power spectrum at large scales was first measured by COBE experiment [l| , and it 
has been confirmed by all subsequent CMB experiments. The spectrum is consistent with a gaussian, and practically 
scale-invariant spectrum, although the latest data from the Wilkinson Microwave Anisotropy Probe (WMAP) Q 
£^ seems to prefer a red-tilted spectrum when no tensor and no other contributions are included, showing a positive 
correlation between the tensor-to-scalar ratio and the spectral index Nevertheless, present data sets at most an 
upper limit on the level of the tensor contribution and non-gaussianity, which even if small, could still be observed by 
the next generation of CMB experiments, like ESA's Planck surveyor satellite 0- 

From the theoretical side, the quest now is for a realistic particle physics model of inflation, and a better understand- 
ing of the inflationary and post-inflationary dynamics. In order to allow for the conversion of the vacuum energy into 
radiation at the end of inflation, the inflaton should couple to other fields. In the standard picture, this process called 
reheating takes place through the perturbative decay of the inflaton into light degrees of freedom, which thermalize 
into radiation. Nevertheless, previous to the perturbative decay, the evolution of the system may be dominated by 
non-perturbative effects as those of preheating [1, i.e., parametric amplification of quantum field fluctuations in 
a background of oscillating fields. Through parametric resonance, field mode amplitudes grow exponentially within 
certain resonance bands in k space, being this a more efficient and faster way of transferring vacuum energy into 
radiation than the standard reheating mechanism. This process and its consequences for the subsequent cosmological 
evolution has been extensively studied in the literature for different kind of models [l(| EH, El, EH ■ The first stages of 
preheating can be studied within linear perturbation theory at first order Q, but soon after the resonance develops 
one would need to improve on the perturbative expansion in order to take into account backreaction and rescat- 
tering effects. Backreaction effects can be partially incorporated by using the Hartree-Fock approximation 0EH- 
However, in order to take fully into account rescattering effects, i.e., mode-to-mode couplings, one has to resort to 
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non-perturbative tools, like lattice calculations [H, [l?], Ell EH, ■ Preheating, although a very fast process compare 
with the cosmological scales, lasting at most only a few e-folds of expansion, may lead to a very ri ch p henomenology 
having to do, among others, with production of massive relics baryogenesis and leptogenesis [U[23[, an d black 
holes production [2(1 HH . 

In addition, in a cosmological framework metric perturbations might also be par ametrically amplified du ring pre- 
heating. This has been studied for example using linear perturbation theory [2^, |2f|, up to second order [2II I28I] . and 
with the Hartree-Fock approximation [29( . A relevant question is whether or not this non-adiabatic excitation of field 
fluctuations leads to the amplification of super-Hubble curvature fluctuations [30l . l3ll [32[ ■ This may happen in the 
presence of an entropy/isocurvature perturbation mode, not suppressed on large scales during inflation, which sources 
the curvature perturbation when parametrically amplified during preheating (33|. Other approach for the study of 
non-linear super- Hubble cosmological perturbations during inflation relies on a gradient expansion, by assuming that 
the spatial derivatives are small compared to the time derivative [U, HH . And more recently it has been developed 
a covariant formalism which would allow to study the full non-linear evolution of cosmological perturbations. It was 
also shown that the non-linear covariant generalization stays constant on super-Hubble scales in the absence of a 
non-adiabatic source [3(1 . 

Besides the amplification of the scalar perturbations, the non-linear nature of the preheating process can also induce 
some level of primordial non-gaussianity in the spectrum, to the level detectable by the Planck mission [2^, [37| ■ In 
addition, it could enhance the tensor perturbations, giving rise to a stochastic background of gravitational waves 
within the reach of the future planned gravitational observatories [38], ]39j] . 

Therefore, preheating after inflation can potentially alter the inflationary predictions. In many cases, like searching 
for non-gaussianity or gravitational waves, we are after very small effects resonantly enhanced during preheating. These 
are by default at least second order in perturbation theory, for which the mixing of the different kinds of perturbations, 
scalar, vector and tensor, is unavoidable. For example, beyond linear perturbation, tensors are seeded not only by 
scalars but also by vectors [10, EH, the latter in turn being generated by the scalars. Aside from preheating effects, it 
has been shown [40] that secondary vector and tensor modes could induced a non- vanishing B-mode polarization in the 
CMB. And this effect becomes comparable or even dominates over that of primary gravitational waves for a primordial 
tensor-to-scalar ratio r < I0~ 6 . This is for example the predicted ratio for models with an inflationary energy scale 
of the order of O(10 15 ) GeV, i.e., one order of magnitude below the typical grand unification scale (42|- Given the 
relevance, for inflationary model building, of the detection of a background of gravitational waves, it is important 
therefore to take into account all possible generating and/or amplifying mechanisms of the B-mode polarization during 
cosmological evolution. In particular, it would be interesting to include in the non-linear simulations of preheating 
the effect of scalar, vector and tensor metric perturbations all together. Although the backreaction effect on the 
fields might be quite small, as expected, they could for example backreact onto each other, affecting the final tensor 
spectrum. 

In most of the previous studies of preheating, fields are evolved in a background metric, i.e., neglecting the effect 
of metric perturbations themselves on the evolution of the fields. It is argued that the backreaction effect of these 
perturbations on the field evolution is small, and the approximation should hold to a good extent when dealing with 
scalar perturbations. Other studies take into account metric variables in a one dimensional system [43| . or reducing the 
system to a one dimensional system by the use of symmetries, either planar [44j or radial [451 ] - These one-dimensional 
studies also focus on scalar metric fluctuations, without vectors or tensors, integrating the system in a lattice with 
N sites. While using lattice techniques allows for a full non-linear treatment of the problem, the procedure has its 
own limitations when dealing with an expanding universe. By discretizing the space and putting the system in a box 
of length L one introduces both a comoving ultraviolet and infrared cut-off in the system. Due to the expansion of 
the Universe, the physical momentum is redshifted and the resonance bands during preheating move towards higher 
values of the comoving wavenumber. Therefore, in order to keep the comoving ultraviolet cut-off larger than the 
effective cut-off for preheating we need a large ratio N/L. On the other hand, if we were interested in super-Hubble 
perturbations, the infrared cut-off would have to be smaller than the Hubble rate of expansion at least at the beginning 
of the evolution, which implies taking L as large as possible. However, limitations on computer memory and CPU 
resources prevend one from working with too big lattices. The advantage of effectively reducing the problem to a 
one-dimensional system lies on the possibility of considering large enough lattices, and therefore super-hubble modes, 
with less computing time cost. We have extended these effectively one-dimensional studies to a 3-dimcnsional system, 
without using any explicit symmetry. We are going to consider preheating of inflaton fluctuations including scalar 
metric perturbations in a cubic 3-dimensional spatial lattice with N 3 sites, and length L. By considering L large 
enough would allow to include super-Hubble modes, but at the expense of quickly loosing the resonance bands during 
the evolution; unless we increment N accordingly, which is not viable in our case. Therefore, as a first step we focus 
on the possible effect of sub-Hubble metric perturbations during preheating. 

In writing down the evolution equations for the scalar field in general relativity, we need first to choose a coordinate 
space-time (t, x) . The choice of coordinates defines a threading of space-time (fixed x) and a slicing into hypersurfaces 
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(fixed t) . In the ADM formalism [46[ , the spacetime line element is given by: 

ds 2 = g^dx^dx" = N 2 dt 2 - llJ {dx i + N i dt)(dx j + N ] dt) , (1) 

with N the lapse function, N 1 the shift vector field, and 7^ the spatial metric. A choice of N and N l fixes the gauge. 
Einstein equations plus the conservation of the stress energy tensor can be split into a set of four constraint equations, 
the Hamiltonian and momentum constraints, and a set of evolution equations for the spatial metric and matter fields. 
The initial conditions for matter and metric variables must satisfy the constraint equations, and cannot be freely 
specified. This is the well known initial- value problem in general relativity [47j] . Once the constraint equations are 
imposed on the initial slide at t — 0, they are fulfilled at any time t by the evolved quantities 1 . We will work in the 
synchronous gauge, which means null shift vector N % = and lapse function N = 1. As a first step, we will only 
follow the evolution of scalar metric perturbations, neglecting vector and tensor modes. That is, we will choose initial 
conditions to ensure that the initial vectors vanishes and the tensors are negligible, and later check to which extent 
they are generated by the non-linear evolution, i.e., to which extent the constraint equations are fulfilled. 

On the matter side, we have studied two different and characteristic types of inflationary models with scalar fields 
minimally coupled to gravity: chaotic models with a quartic potential A<& 4 , and standard supersymmetric hybrid 
model. Preheating in both of them has been extensively studied in the literature, including lattice simulations. For 
the chaotic quartic model the resonance starts in a single, narrow resonance band [10| at a fixed value of the comoving 
wavenumber, but later on due to rescattering the resonance spreads to several bands in momentum space, until finally 
the spectrum is smoothed out [l6[ . In hybrid models due to the tachyonic instability at the end of inflation [l2j the 
lower momentum modes are quickly amplified, and in a few oscillations of the background fields the whole process 
ends 0,0. 

In section II we give the set of equations for metric and field variables to be integrated in the lattice, and discuss 
the issue of the initial conditions. In section III we present the results for the chaotic model with a quartic potential; 
results for the hybrid model are given in section IV. In both cases, we compare the evolution with and without scalar 
metric perturbations. The summary of our results is given in section V. 



II. EVOLUTION EQUATIONS AND INITIAL CONDITIONS 

Our choice of gauge reduces the metric given in Eq. (Q]) to the line element: 

ds 2 = e 2a{t ^dt 2 -e 2 ^%dx l dx\ (2) 

where we have defined the lapse function N = e a , and the spatial metric 7^ = e 2/3 jij, with Det(ff) = 1. Although 
we are going to work in the synchronous gauge and later set N = e a = 1, we will keep the dependence on the lapse 
function explicit in this section for the sake of generality. The matter content is given in general by a set of scalar 
fields x), with the Lagrangian and stress-energy tensor: 

£ = iv^'v*! - V(*j) , (3) 

V - Vn&VuQi - \g^ (V CT $ J V CT $i + 2V($/)) , (4) 

where the subindex "I" counts the no. of fields and summation is understood over repeated index; V($/) is the 
potential for the fields considered, and V M the covariant derivative for the metric g^. With this choice, the evolution 
equation for each scalar is given by: 

+ 3f3'^'j - + Vi - D t aD l <S>! = , (5) 

where now Di is the spatial covariant derivative given by the spatial metric 7^ , and we have defined: 

^=e- a d t ^j I3' = e~ a d t (3. (6) 

In the covariant formalism [49j], these are the covariant time derivative along the unit four velocity = dx^/ds = 
(e _Q , 0,0,0), with /' = u^V^f for any scalar quantity. The volume expansion is then given by V^m^ = 3/3', and the 



1 Although the discretization procedure will invariably introduce deviations from the constraints during the evolution [4§| . 
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acceleration vector by a 11 = u v V „it M = (0, D l a). Given the four vector u^, the stress-energy tensor can be decomposed 
as: 

T^y = (p + P)Ufj,u u + g^P + q^Uu + q v u^ L + , (7) 

where p, P, and tt^ are respectively the energy density, pressure, momentum and anisotropic stress in the frame 
defined by u M , and given by: 



p = 




p = 


I 


q% = 




n« = 





\4>i) , (8) 

-vm, (9) 

(10) 

(11) 

Einstein equations give the evolution of the spatial metric components, plus 2 constraint equations. With our choice 
of gauge, the evolution equation for (3 reduces to: 



3/3" + 3(/3') 2 - Di&a = - - (p + 3P) - A} Af , (12) 



where k = 8itGn = 1/™%, with mp = 2.42 x 10 18 GeV the reduced Planck mass, and Aj is the traceless part of the 
extrinsic curvature Kj A6], which in our case is given by 

K i = tl^ dtlkj= psi + A i j , (13) 

where Aj = 7 7j.-/2. The latter contains the vector, tensors and the traceless scalar mode, while /?' gives the local 
expansion rate. The equation for Aj reads: 

A! + 2,(3' A) = e- a ^D l Dje a - ^D k D k e a J ~ (^ 3) R) - ^ (3) i?^J + K (ll) - ^IlJ , (14) 

where ^R), ^R are the Ricci tensor and scalar with respect to the spatial metric. For example for a flat spatial 
metric with an inhomogeneous scale factor, 7^ = e 2/3 <5y, we have: 

{3) R) - | (3) i? = e- 2!i {tfpdjp - d%f3 - |(d fe /3d fe /3 - d k d k ^ , (15) 
while the field source term is given by: 

nj - = e - 2f3 ^8 i &d j & - ^-a fc $9 fc $^ . (16) 

The Hamiltonian and the momentum constraint are given respectively by: 

60'" -^R-A)A{ = 2k P , (17) 
2Di0-DjA\ = -K&jDi® 1 . (18) 

The full set of equations ([5]), (fl^j) and (fT4"|) have been integrated in a lattice in Ref. [5(| for a chaotic inflation 
model, where it was shown the viability of inflation starting with inhomogeneous initial conditions. We will instead 
follow the evolution starting just at the end of inflation. However, due to the complexity of the problem, as a first 
step we will only keep the effect of the scalar metric variable j3, assuming that vector and tensors are subdominant 
and negligible. Afterward we will cross-check the consistency of this approximation. Therefore, keeping only the (3 
dependent terms the final set of equations to be integrated is given by: 

+ 3l3'& I -e- 20 d l d l ® I + V I -e- 2 ' 3 (d l a + d l f3)d l <P I ~ 0, (19) 
3/3" + 3(/3') 2 - e-^i&'diO. + d'ad.a + d^d.a) ~ ~{p+ 3P) . (20) 
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The initial conditions for <!>, /3, 0' must be compatible with the constraints. For the fields we follow the standard 
procedure, with the quantum field theory being replaced by an equivalent classical field theory. Equivalent in the 
sense that the expectation values of quantum variables are equal to the average values of the classical ones. The 
classical to quantum transition takes place during inflation, when interactions of the quantum fields can be neglected. 
Therefore, the field initial conditions for the subsequent classical evolution can be obtained from the linear theory of 
quantum-to-classical transition [5ll | . 

The field is expanded in Fourier modes in a spatial lattice of volume L 3 with periodic boundary conditions. The 
zero mode corresponds to the homogeneous value of the field given by its average value in the lattice, Q(t) = ($(t, a;)). 
The non-zero modes are expanded in fc-space as: 



= i fc)e*% (21) 

k 

dt*i(t,x) = ±J2 d Mt,k)e ikx . (22) 

k 

The initial vacuum state q(0, k) corresponds to a complex Gaussian distribution with a random phase for each mode 
k, and root mean squared 

\q(0,k)\ rms ~ — , (23) 
v 2^^ 

where uuk — Vk 2 + Vn is the initial frequency of each k-mode. The initial time derivative for the field is then given 
bydtq(0,k) = -w k (0jq(0,k). 

Having set the initial field's values, we have to deal now with the initial conditions for the metric variable, such 
that the constraint equations are fulfilled at t = 0. Those are given by: 

3/3' 2 -e- 2f3 {d k pd k l3 + 2d k d k l3) = K p, (24) 

2di(3' = -K&jdi® 1 . (25) 

The main restriction comes from the momentum constraint: without including metric vector variables, the RHS of 
Eq. (|25p must reduce to a gradient term at least at t = 0. However, this is not the case in the synchronous gauge 
when studying single field models with initial conditions for the field as given in Eqs. (f2Tj) and (|22|) . In order to be 
consistent with our assumption of negligible vector and tensor modes, we can take instead an homogeneous initial 
profile for the field velocity, given by its background value, i.e, <£>' = (dt$(0,x)). This will set our choice of initial 
conditions for the chaotic model studied in section III. Another option often used in the literature could be to add 
an extra, non-interacting scalar field [44| . Its initial velocity can be adjusted such that the momentum constraint 
is trivially satisfied [52j . The energy density of the non-interacting extra field would behave as radiation due to the 
gradient contribution in Eq. l[8]). and setting it to be negligible at t = 0, it would remain so during the evolution. 
We will use instead this choice for the hybrid model in section IV, where the role of the "extra" field is played by the 
waterfall field already present in the model. 

Having adjusted the initial field velocities in order to fulfill the momentum constraint, Eqs. and ([25]) provide 
the initial values of f} and /3' at t = 0. For one-dimensional systems (or effectively one-dimensional) this procedure 
is consistent with neglecting vector and tensor perturbations, and the constraints are preserved by the evolution. 
However, in general this is not the case in 3 dimensions without additional symmetries, and vectors and tensors will 
be induced during the non-linear evolution. Therefore, the departure from zero of the momentum constraint during 
the evolution may give us an estimation of the amplitude for example of the vector perturbations generated during 
preheating 2 . On the other hand, the scalar perturbations we consider will contribute to the source term for tensor 
perturbations, i.e, the RHS of Eq. (|14p . and we can compare this to the pure anisotropic scalar field source term, as 
those considered in (39j . and check whether or not they may become comparable during preheating. 



III. CHAOTIC MODEL 



We first study the chaotic inflation model with potential V = A$ 4 /4. In this model, inflation takes place for 
values of the background field larger than mp, and ends approximately when <j)(t) = ~ 2y/3mp. Following 



2 Tensor perturbations, being transverse and tracclcss, do not contribute to the momentum constraint. 
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FIG. 1: Left hand plot (a): Spectrum of irmaton field fluctuations, normalized by the initial value of the field, a(t) 2 |<&fc| 2 , at 
different times, (a) t = 0, (b) t = 10 3 , (c) t = f0 4 , (d) t=5 x f0 4 , (e) t = 10 s , (f) t = 5 x f0 5 . Solid lines are the spectrum 
with scalar metric fluctuations included. Dashed lines show the spectrum without metric fluctuations. Right hand plot (b): 
Spectrum of scalar metric fluctuations a(t) 2 \{3k\ 2 at different times. We have defined a(t) 2 = (e 13 ) 2 . 



Ref. [lq , we start evolving the field when oscillations begin, at a slightly smaller value 4>q = (^(O)) = \/3mp. 
Approximately this is the value at which the conformal time derivative of the field vanishes. We rescale the field by 
its initial value (j)o, and time will be given in units of l/\/A</>o. In program units, the initial background field velocity 
is (j)Q = ($(0)) = — 1/V2- The comoving wavenumber is given in units of \/A</>o. With this rescaling the coupling A 
does not appear in the equations of motion, but it rescales the initial value of the fluctuations (the initial spectrum). 
Given that the development of the resonance for the inflaton fluctuations is not very sensitive to its value, we will take 
A = 1CP 4 for convenience, although in a realistic model it should be A ~ O(10 -14 ) in order to match the amplitude 
of the primordial spectrum with the observed value. 

The initial values for the field fluctuations are given by Eqs. (j2"Tj) . and (|2"5)l . and we adjust the initial value of the 
field velocity in order to fulfill the momentum constraint at t = Q, such that 

$'(0) = (dMO)) . (26) 

Finally, the initial profiles for the local scale factor e 13 and expansion rate /3' are obtained from Eqs. (pM)) and (|25p . 
with (e^ ') = 1. 

We have run the simulations in a lattice with TV = 64 and L = 107T, and periodic boundary conditions for the field. 
In program units, the smaller comoving wavenumber is fc = 0.2, and the initial value of the average expansion rate 
H = (/3'(0)) ~ 1/V2. That means that at the beginning of the simulation there are a few field modes that have 
crossed the horizon before inflation ends, with k < a^Ho, with ao — (e^ )), and have not yet reentered. For those 
modes the amplitude of the spectrum has frozen around the time k = a*iJ», and therefore the initial value would be 
slightly larger than Eq. (|2l?|) by a factor a /a* [53[ , given by 

kfo 3/«o fc _\ 
a* i \a* / 

with k given in program units. In Fig. |T]a) is shown the spectrum of field fluctuations |a(i)$fe| 2 at different 
times, where a(t) — (e' 3 ) is the averaged scale factor. We compare the results obtained when including scalar metric 
fluctuations (solid lines) with those obtained by evolving the field with a background metric (dashed lines). At low 
values of the comoving wavenumber k it can be seen the initial enhancement of the amplitude for the superhorizon 
modes, which is further amplified by the parametric resonance. But the main features of the resonance in this model 
are unchanged by the inclusion of superhorizon modes and/or metric fluctuations, and have been well established in 
the literature [lO, EB, E3, [H3|- First, field fluctuations grow in a narrow resonance regime, with the resonance peak 
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aCt) a(t) 



FIG. 2: (a) Left hand plot: average total energy density (p) versus the average scale factor a(t) (black solid line). Also 
shown are the average gradient contribution of the field Gradient={e -2 ^| V$| 2 /2), and that of the scalar metric fluctuation 
Grad(/3) = (e~ 2/3 \Vf3\ 2 /2). The dashed grey (red) line shows the gradient of the field without metric perturbations, (b) Right 
hand plot: Variance of the field (solid bottom line) and scalar metric fluctuation (solid top line). The field variance without 
metric fluctuations (dashed line) is also shown for comparison. 



located at k ~ 1.27. At later times, rescattering effects lead to the appearance of multiple peaks, until finally the 
spectrum becomes smooth at small momentum, with an ultraviolet cut-off increasing in time. It would be followed by 
a long stage of thermalization, during which the transfer of power from low to high momentum continues in the free 



turbulence regime [54j. However, this stage is difficult to study numerically because of the limitations of the lattice 
cut-off, fixed by the choice of the lattice size, which soon becomes smaller than the physical cut-off which is increasing 
with the scale factor. 

The main difference when including metric fluctuations can be seen at smaller momentum during the first stages of 
the resonance. The initial superhorizon modes tend to be amplified at the same time than the main resonance peak, 
and this effect is enhanced by the metric fluctuations. This was already observed in Ref. [44| in the 1+1 dimensional 
system. Our simulations show that the effect prevails in 3+1 dimensions. Nevertheless, at later times the tendency 
seems to be reversed, and the amplitude of the long wavelengths modes starts to be smaller than when integrating 
the system with a background metric. This might be an indication that due to the initial increase, the transfer of 
power from long to short wavelengths starts before and it reaches sooner the stage of free turbulence. 

The spectrum of the scalar metric fluctuation \f3k\ 2 is shown in Fig. |T]b). The initial spectrum (top solid line) 
behaves like 

|/3 fe (0)| 2 «^f|$ fc | 2 , (28) 

for k > 1, showing the same enhancement at low k (superhorizon) than that of the field. With time, the spectrum is 
redshifted like the inverse squared of the average scale factor a(t), but it retains the resonance peak structure of the 
field. The final spectrum (dot-dash-dashed line) is given by the initial spectrum redshifted by a factor a(t)~ 2 , but 
with the initial power at low k smoothed out, and showing the same cut-off at large k. The field amplitude in this 
model is redshifted like a(t) , but the rescaled field mode a(£)$fc exponential increase due to parametric resonance. 
This is not the case for the metric fluctuation, where the fluctuation is not really parametrically amplified. 

When looking at the averaged values of fields, energies and variance, the inclusion of scalar metric fluctuations has 
a negligible effect. In Fig. (0a) is shown the average total energy density (p) with respect to the average scale factor. 
The system behaves as radiation with (p) ex a(t)~ 4 , and the average value is the same than that obtained without 
metric fluctuations. There is at the beginning a slight difference in the gradient contribution of the field (e~ 2 ^|V<I>| 2 ), 
but it practically disappears once the resonance broaden due to rescattering effects 3 at around a(t) ~ 50. The gradient 



For comparison with the times at which the power spectrum is plotted, we have the standard time dependence of a radiation dominated 
universe, with a(i) ~ y/i. 
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FIG. 3: Average momentum constraint Cm (see text), with (solid line) and without (dashed line) metric fluctuation. 

term of the metric fluctuation (e — 2/3 1 V/3| 2 ) initially follows that of the field, redshifted as a~ 4 , but the resonance in 
the field is not strong enough to keep this tendency and soon it starts decreasing faster as a~ 6 . Side by in Fig. ((2]b) 
we have included the variance of the field (<I> 2 ) — ($) 2 (solid and dashed lines) and metric fluctuation (0 1 ) — ((3) 2 (solid 
line). Although the variance is initially larger with an inhomogeneous scale factor /?, the difference again practically 
disappears when the resonance develops. The dispersion in j3 is larger at the beginning than that of the field, but 
again it is less affected by rescattering effects. 

Finally, we have checked to which extent metric vectors fluctuations can be neglected in the field and scalar metric 
evolution, by checking the momentum constraint Eq. (|25|) . Although initially we set the vectors and tensor fluctuations 
to zero, they are sourced by scalar perturbations. The failure to satisfy this constraint during the evolution can be 
interpreted as an indication of their presence. This can be seen in Fig. ([3]), where we have plotted, as a function of 
the scale factor, Cm defined by: 

Cm = W + ^'OM 2 ) 1 ' 2 , (29) 
which shows to which extend Eq. (|25|) is verified. This will accounts for the terms neglected in Eq. ([18]) . such that 

C m = {\\d 3 A{\Y' 2 , (30) 

We include for comparison also the result without scalar metric fluctuations, i.e., with an homogeneous scale factor. 
Although when including (3 the constraint Cm seems to be better fulfilled by and order of magnitude with less 
dispersion, in either case it just follows the behavior of the field fluctuation, with Cm decreasing initially like (<fr'di&) oc 
a(t)~~ 3 until the exponential increase due to the resonance of the field, fndeed the behavior of Cm indicates that also 
the gradient of the local Hubble parameter follows that of the field, with (di/3') oc a(t)~ 3 instead of being redshifted 
like a(t)~ 2 . We have checked that the behavior of Cm is not a numerical artifact due to discretization by integrating 
the system for other choices of L and N . 

Therefore, vectors (and by extension tensors) could be generated with an initial amplitude similar to that of scalar 
fluctuations. However, due the behavior of the gradient term for (3 in Fig. ([Ha), the source term in their evolution 
equation (TH|) will be soon dominated by the contribution of the scalar field. Given that the resonance is not strong 
enough to be fully transferred from the fields to the scalar metric fluctuations, we might anticipate that the same 
would happen for the other kinds of perturbations. Nevertheless, scalar metric fluctuations affect the field spectrum 
at the start, although their effects are quickly redshifted, and vectors and tensors may lead to a similar effect. It would 
be interesting to study to which extent metric fluctuations apart from the scalar ones help to enhance the power of 
the lower modes, whether this process last for longer or just speed up the transfer of power through the spectrum, 
but this is beyond the scope of the present work. 
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IV. HYBRID MODEL 



Preheating effects strongly depend on the model considered, and models with an stronger resonance will lead to 
different conclusion also about metric perturbations. Because of that, we study now a different kind of model, a 
standard suspersymmetric hybrid model of inflation, where it is known that the field resonance is stronger and driven 
in what is called the tachyonic regime [H, EE Ell • The potential is given by: 

V = V + ^N 4 + 5 2 ($ - $+)($ - <P-)N 2 + X - m \<$> 2 , (31) 

where as before $ is the inflaton field, with a small mass rr£, and N is the waterfall field which triggers the phase 
transition at the end of inflation. Inflation takes place in the false vacuum, with (N) ~ 0, and V ~ Vq. When 
the inflaton field overcomes the critical value the squared field dependent mass for N becomes negative, due 
to the tachyonic instability field fluctuations grow exponentially, and both fields move towards the global minimum 
of the potential at O = ($£ + $^)/2 and N = (<f>+ - <f>~)/\/2. The potential in d3TJ) reduces to the standard 
supersymmetric hybrid model when $~ = — but it also includes other supersymmetric hybrid potentials when 
the soft trilinear term is included [Hj]. Hereon we will take = 2</>o, and <5>~ = 0; other choices of these parameters 
will only shift the background values of the fields at the global minimum, without affecting the evolution of fields 
and field fluctuations. The vacuum energy Vq driving inflation is adjusted such that the potential energy vanished 
at the global minimum, i.e., Vo ~ g 2 N^/A = g 2 <pQ. As before, in order to eliminate the coupling from the evolution 
equations we rescale the fields by </>o, the potential and energy density by g 2 4>Q, time will be given in units of l/(g<f>o), 
and wavenumber in units of g4>a- The value of the coupling then sets the initial value of the fluctuations. 

Working only with the background fields during inflation, i.e. <f> = (<£>), the COBE amplitude of the primordial 
spectrum is given by: 

(32) 



2ttJ 47rV3r/0 \m P 

with 770 = rr&/ (3H 2 ) one of the slow-roll parameters, which also gives the tilt of the spectrum with the spectral index 
ns — 1 + 2?70. Depending on the values of the parameters in the potential, the coupling g and the mass m^, the 
scale of hybrid inflation can range from the unification scale, with ^fg~4>o — O(10 15 GeV), down to the supersymmetric 
breaking scale 0(1 GeV). Numerically, a too low value of the scale means a lower value of the Hubble rate of expansion 
and requires more integration time until the fields start oscillating. Therefore we will work here with models near the 
GUT scale, and start with O = 0.005m P ~ 1.2 x 10 16 GeV, g = 0.01, ^ = 0.05 and V Q 1/4 ~ 1.2 x 10 15 GeV. With 

these values of the parameters 4 we have P^/ 2 ~ 5 x 10~ 5 , and a very small tensor-to-scalar ratio r ~ 8 x 10~ 6 . 

We begin the evolution of the system just before the end of inflation, i.e., some fraction of e-fold N e ~ 0.05 before 
the end with the value of the background inflaton field just above the critical value ($(0)) = $+e n < t>N ' i . We use the 
slow-roll condition to set the initial value of the background velocity given by ( < i ) (0)) ~ —2r\^(f>Q / 'mp in program units. 
The background value of the waterfall field (N) would be around zero, although fluctuations very quickly drives this 
to a non zero value. Numerically, it does not make much difference if we take a non vanishing initial value, and we 
have set (N(0)) ~ 10 -8 and (iV(0)) = 0. The initial values of the inflaton fluctuations are given by Eqs. (|2Tp. (j2"2"|) and 
(f2"3"| . But now, in order to fulfill the momentum constraint at t — 0, we adjust the initial value of the N fluctuations, 
identifying this field as the extra field considered in [H^]. Therefore, we have taken: 

N(0) = HO) - ($(0)) + (N(0)) (33) 
AT(0) = ($(0)) - 4(0) (34) 

and the momentum constraint at t = is given by: 

28,13' = -K($(0))8i$(0) ■ (35) 

With the choice of parameters given above, the Hubble rate of expansion is rather small and practically the expansion 
effects are going to be negligible during the tachyonic resonance. In particular, in program units we have -ff(O) = 



4 Strictly speaking we also have a blue tilted spectrum with a rather large tilt, n$ — 1.1, excluded by observations 5, 6] when the tensor 
contribution is negligible. We could always take a smaller value for both the coupling and r)^ and satisfy all constraints. Our choice of 
?70 for this value of the inflationary scale is motivated by numerical considerations, in order to shorten the first stage of the integration 
from the end of inflation up to the point when fluctuations starts growing. 
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FIG. 4: (a) Left hand plot: average total energy density (p), potential (V), kinetic (i> 2 + N 2 ) /2, and gradient (e 2,3 ((V$) 2 + 
ViV) 2 ))}/2 energy density, versus time (solid lines). The average gradient contribution of the scalar metric fluctuation 
<e- 2/3 |V/3| 2 )/2 (solid grey/ green line), scale by a factor of 100, is also shown. The dashed lines show the contributions without 
metric perturbations, total energy density, potential and field gradient energy density. Field quantities are given in units of 
the initial energy density g 2 (f>o, that of the metric fluctuation in units of g 2 (f>o/irip. (b) Right hand plot: Variance of the fields 
(solid lines), and scalar metric fluctuation scaled by a factor of 100 (solid grey/green bottom line). The field variance without 
metric fluctuations (dashed lines) is also shown for comparison. The dot-dashed lines are the average field oscillations scaled 
by a factor of 10 -2 . Field quantities are given in units of <f>%, that of the metric variable in units of {cj>o / mp) 2 . We have taken: 
0o = 0.005 m P , k = 0.01, jfa = 0.05. 

(/3'(0)) = (j> /(V3m P ) ~ 3 x 10~ 3 . This also means that at the beginning of the evolution all modes considered are 
well inside the Hubble radius. 

As before, we integrate the system in a cubic lattice with N 3 = 64 3 sites and L = 10tt, and periodic boundary 
conditions for the fields. Preheating in this model proceeds through the exponential growth of the waterfall field 
fluctuations due to the tachyonic instability in the potential, which makes those of the inflaton to grow at the same 
rate [T3.[l8j. and those of any other field coupled to it [23, 55]. This lasts only a few oscillations of the fields, depending 
on parameter values, during which the energy density ends more or less equally distributed among potential, kinetic 
and the gradient contribution of the fields. The tachyonic instability ends when the field fluctuations of $ and N 
render the effective squared mass for N always positive during the oscillation, i. e, when the variances of the fields 
become O(tf>o)- Following this there would be a period of bubble collisions [19j . and then a stage of turbulence with 
the transfer of momentum towards the ultra-violet, and finally the thermalization of the system. As with the previous 
chaotic model, we concentrate here in the first two stages without really following the system into the turbulence 
regime. These can be seen in Figs. (@|, where in the LHS we have plotted the total energy density, potential, kinetic 
and gradient, and in the RHS the variance of the fields. All quantities are given in program units, i.e., energy densities 
are normalized to the initial vacuum energy Vb, and the variances to the field value 4>q. Wc include for comparison 
the results with no metric perturbations. There is no qualitative difference between both cases, except for a slight 
shift in time. In the RHS plot we have also included the averaged gradient term for the scalar metric perturbation 
(e _2/3 |V/3| 2 ), rescaled by a factor of 100. This term is always negligible with respect to the field contributions, and 
therefore metric perturbations has little or no effect on the evolution of the fields. However, in this case the metric 
fluctuations are also exponentially amplified, following the same resonance pattern than that of the fields. This can 
be seen in Fig. ([5]). On the LHS it is shown the evolution with time of the waterfall field spectrum | TV^ | 2 , and in 
the RHS that of |/3fc| 2 . Both spectra behave during the resonance as e 2MfcAt , with a growth index fik ~ 0.3 for the 
lowest modes, slightly larger for /3k- In addition, the size of our lattice L = lOir is small enough to keep the physical 
ultraviolet cut-off in momentum for the fields, with the largest momentum mode kept outside the resonance band and 
being hardly enhanced. But this is not the case for the scalar metric fluctuations, for which all modes are preheated. 
The initial scalar metric spectrum is derived from the energy constraint Eq. (|24p . and therefore |/3fc| 2 oc O(((f>o/mp) 2 ). 
The exponential growth is not enough, being only slightly larger than that of the field, to make the scalar metric 
contribution in the evolution equations comparable to that of the fields. In particular in the evolution equation for 
the traceless metric components, Eq. (1141) . the source term due to the scalar metric fluctuations would be of the order 
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FIG. 5: (a) Left hand plot: Time variation of the waterfall field spectrum |iVfc| 2 , in program units, (b) Right hand plot: Same 
for the scalar metric fluctuation |/3fe| 2 . We have taken: (j>o = 0.005mp, k — 0.01, 77^ = 0.05. 

0(e~ 2/3 |V/3| 2 ), and therefore suppressed by a factor O(((f)o/mp) 2 ) with respect to the field source term. For the same 
reason, the traceless components A\ would be of the order Ke~ 2/3 |V<I>| 2 ~ 0((4>o jrap) 2 , and their contribution to the 
evolution of fields (O, and scalar metric fluctuation p2[) . negligible. 

However, given that the main suppression of the scalar metric terms with respect to those of the fields is due 
to factors of the order of O((f>o/mp), we have checked whether this scenario changes when increasing the rate of 
expansion, and/or decreasing the value of and the coupling g. In particular, we want to check whether or not the 
source term due to metric perturbations in the evolution equation for the traceless modes Eq. (|14|) can be comparable 
to that of the fields. In Fig. ([6]) we have plotted the averaged value of the field's gradient and that of the scalar metric 
variable /3, normalized by the initial vacuum energy density, for different choices of the parameters. In all of them 
we have increased the rate of expansion by and order of magnitude, with tfio/mp — 0.05 and then Hq ~ 3 x 10~ 2 in 
program units. Keeping first the other two parameters the same than in the previous plots, = 0.05 and g = 0.01, 
we have already that the scalar metric contribution becomes comparable to that of the fields (solid lines). In this case, 
due to the larger Hubble expansion rate, tachyonic preheating ends during the first oscillation of the fields. Decreasing 
the value of = 0.01, i.e., making the inflaton potential flatter, slightly delays the beginning of the oscillations and 
the start of the resonance (dashed lines), but qualitatively the behavior of the gradient terms remain the same. Same 
when decreasing the value of the coupling, with g = 0.001 and j],p = 0.01 (dot-dashed lines). By lowering the value 
of the coupling we decrease the initial value (relative to the total energy density) of the fluctuations. Until the fields 
start oscillating, they behave approximately like massless fields, and the gradient terms are redshifted like a(t)~ 4 , 
with the effective scale factor a(t) = (e^). When further decreasing 77^ (dot-dot-dashed lines), the initial value for the 
gradient terms is further decreased, but the relative enhancement during the resonance remains the same. However in 
this case the gradient of the fields decreases faster after the tachyonic resonance than that of /3, such that the latter 
soon dominates after the resonance. 

Nevertheless, when increasing the scale of inflation we loose the physical cut-off in our lattice, in the sense that 
the highest momentum mode is further enhanced once tachyonic reheating ends. In this case, it just take less than 
one oscillation for the variance of the fields to grow enough to render the effective squared mass of the waterfall field 
positive. From the numerical simulations we have that for g = 0.01 and 77^ = 0.05 this happens when ait) ~ 2.3; for 
g = 0.01 and 77^ = 0.01 when a(t) ~ 3.5; and for g — 0.001 and 77^ = 0.01 when a(t) ~ 5.6. Therefore, practically all 
modes have been redshifted inside the resonance band by the time the fields hit first the global minimum. In order 
to check the dependence of the results on the lattice cut-off, we have run a simulation with N = 128 and L = lOw, 
i.e., doubling the maximum value of the comoving momentum, when g = 0.001 and 77^ = 0.01, such that still at the 
end of the resonance there are modes left outside the resonance band. This is shown in Fig. ^ by the dotted lines. 
The values of the gradients behave the same than before, and they do not differ by more than a factor 0(2). 

In any case, during the short period of the resonance all long wavelenghts modes arc excited, and immediately after 
it seems to start a very effective transfer of momentum from the infrared to the ultraviolet. This can be seen in Fig. 
([7]), where as an example we have plotted the spectrum of the waterfall field liV^ 2 (LHS plot) and that of \(3k\ 2 (RHS 
plot), for 0o — 0.05m p , g — 0.01, 770 = 0.01. On the latter, it can be seen more clearly the point at which the tachyonic 
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FIG. 6: Comparison of the gradient of the fields (black lines), Gradient=(e~ 2/3 ((V$) 2 + ViV) 2 )))/2, and that of the scalar 
metric fluctuation (grey/red lines) Grad(/3) = (e~ 2 (V/3) 2 )/2, normalized by the initial energy density, for different parameter 
values and N = 64, L = 107r: g = 0.01 and = 0.05 (solid lines); g — 0.01 and r/^ = 0.01 (dashed lines); g — 0.001 and 
?70 = 0.01 (dot-dashed lines); g — 0.01 and 77^ = 0.001 (dot-dot-dashed lines); in a lattice with N — 64 and L — 10tt. We also 
include for comparison the results obtained in a lattice with N = 128, L = IOtt, and g = 0.001, = 0.01 (dotted lines). For 
all we have taken cj>o = 0.05m p . 

resonance ends at around t ~ 30 when the highest modes stop to be exponentially amplified, and the restart of the 
enhancement of the highest modes at around t ~ 40. Same happens for the other model parameters considered in Fig. 
((SJ), at different time intervals. Given that the growth of the gradients practically ends with the tachyonic resonance 
as expected, we do not expect their evolution to be much affected by the lack of high momentum modes. As in the 
standard case without metric fluctuations, the resonance will be followed by the turbulence regime, which study we 
do not pursue here but it would be necessary in order to see the scaling followed by the metric perturbation terms 
there. In Fig. © the j3 terms decay as fast as that of the field, if not slower, but this tendency remains to be further 
checked. But we would like to stress that, even if non-linear metric perturbations are not preheated enough to affect 
the evolution of the fields, they may become comparable to the field source term in their own evolution equation, 
depending mainly on the scale of inflation. Our first study here indicates that this is the case when <fio/mp > 0.05, 
which sets the effective rate of expansion during the numerical evolution, and it does not depend much on the value 
of the coupling <?, neither on that of the scale of inflation. For example a hybrid model with (j>o/rnp = 0.05 and 
770 ~ 0.01, we will get the right amplitude of the primordial spectrum with g ~ 10~ 4 , which means a scale of inflation 

Vo /4 ~ 10 15 GeV. The same scale is obtained with g = 0.01 and <po/mp = 0.005, but in the former model we may 
expect preheating of the metric perturbations to be non-negligible, whereas in the latter we have checked that this is 
not the case. 

V. SUMMARY 

We have studied the first stages of preheating following inflation in a 3 dimensional spatial lattice, including both 
field and metric perturbations. This extends previous one dimensional studies of the problem as those of Refs. 
[43I [44], [4||. We have worked in the synchronous gauge, starting the evolution of the system immediately after 
inflation. The set of evolution equations for fields and metric variables are given by Eqs. (J5J), (|12[) and l|14p. During 
inflation, vector modes decay, and tensor modes are subdominant with respect to the scalar ones. Thus we have 
taken them, vectors and tensors, initially to vanish. This means that we have to take initial profiles for the field and 
field velocity consistent with the momentum constraint Eq. (|18[) . The non-linear evolution of the system couples all 
kind of modes to each other, and any of them will be inevitably produced. Nevertheless, as a first step we have only 
followed the evolution of the scalar metric perturbations and their effect on the field fluctuations. We have studied two 
different kind of models of inflation, a chaotic model of inflation and a standard hybrid one, which are representative 
of different patterns for the parametric resonance. In the former, the resonance starts in the narrow resonance regime 
with field mode fluctuations produced in a narrow range of comoving wavenumber Q . Later on, rescattering effects 
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FIG. 7: (a) Left hand plot: Time variation of the waterfall field spectrum |iVfc| 2 , in program units, for <po = 0.05mp, g = 0.01 
and r/$, = 0.01. N = 64 and L = Wit (b) Right hand plot: Same for the scalar metric fluctuation \Pk\ 2 . 



re-distribute the resonance among several bands, for higher and lower comoving momentum values [HI, [l?]]- In the 
hybrid model the resonance is driven by the tachyonic instability in the waterfall field, and due to that it results in a 
more explosive and faster production of particles [121, [l8|, [l9| . In both cases, scalar metric perturbations hardly affect 
the evolution of field fluctuations, and the evolution of the field resonance is practically unchanged. For the chaotic 
model of Section II, however, there is an initial enhancement of the field power at lower wavenumbers, that was also 



observed in Ref. [441 ) . In addition, we have also included the effect of superhorizon modes in our simulations, with 



larger initial amplitudes. Looking at the end of the integration range in time, there might be also some indication 
that this initial transfer of power towards lower modes helps in entering earlier the stage of free turbulence, where the 
effect is reverse. But given the size of our numerical simulations we cannot follow the system fully into that regime, 
and further studies would be required in order to confirm this tendency. 

Although the field resonance does not feel the presence of metric fluctuations, even in the non-linear regime, the 
scalar metric perturbations follow the same resonance pattern than the field, as seen in Fig. (jTJb) for the chaotic 
model, and Fig. (0b) for the hybrid model. In the case of the chaotic model the amplitude of metric perturbations is 
not really amplified above the initial value, and the power spectrum of the metric fluctuations simply reproduced the 
peak structure seen for the field. This means that for the chaotic model metric fluctuations are not really enhanced 
by the parametric resonance, they do not become comparable to the field source terms in the evolution equations, 
and we can safely neglect them. On the contrary, in the hybrid model the scalar metric perturbations clearly feel the 
tachyonic resonance of the waterfall field, in a similar way to any other field coupled to it like the inflaton field. The 
scalar metric initial amplitude is amplified by approximately the same factor ~ e fJ ' kAt than the amplitude of the field 
fluctuation, with ~ 0.3 being the growth index. Although metric fluctuations never become large enough to affect 
the field evolution, they may give rise to a source term comparable to that of the field in the evolution equation for 
the traceless modes, Eq. (I14p . depending on the model parameters g and (fro/nip, with Vq^ 4 = ^fg4>o- For example, 

this is the case for g = 10~ 4 and 4>a/mp — 0.05, i.e., for a scale of inflation ~ 10 15 GeV. Whether this term helps 
to enhance or not for example the production of gravitational waves during preheating [H, [3t| , requires the study of 
the full set of evolution equations for the metric perturbations, which is beyond the scope of this paper and is left for 
further study. 
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